Thermostatistics of overdamped motion of interacting particles 
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We show through a nonlinear Fokker-Planck formalism, and confirm by molecular dynamics sim- 
ulations, that the overdamped motion of interacting particles at T = 0, where T is the temperature 
of a thermal bath connected to the system, can be directly associated with Tsallis thermostatistics. 
For sufficiently high values of T, the distribution of particles becomes Gaussian, so that the classical 
Boltzmann-Gibbs behavior is recovered. For intermediate temperatures of the thermal bath, the 
system displays a mixed behaviour that follows a novel type of thermostatistics, where the entropy 
is given by a linear combination of Tsallis and Boltzmann-Gibbs entropies. 
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Nonlinear Fokker-Planck equations (FPE's) [l| are fre- 
quently employed to represent macroscopically physical 
and chemical systems displaying anomalous diffusion be- 
havior 0. Scientifically and technologically important 
examples of such systems include, among others, the flow 
through porous media 0] , the dynamics of surface growth 

the diffusion of polymer- like breakable micelles 
and the dynamics of interacting vortexes in disordered 
superconductors A typical nonlinear FPE may be 

written in the general form [l], Q , 



dP 
~dt 



d[A(x)V(P)] 
dx 



d_ 

dx 



dP 

dx 



(1) 



where the external force A(x) is associated with a po- 
tential 4>{x) [A(x) = —d<f>(x)/dx], and the analyticity of 
the potential as well as the integrability of the force are 
assumed to hold in all space. The functionals ^[P(x,t)] 
and Sl[P(x,t)] satisfy requirements of positiveness, inte- 
grability, and differentiability with respect to P{x,t) @. 
Moreover, in order to preserve the probability normaliza- 
tion for all times, one should impose the probability dis- 
tribution, together with its first derivative, as well as the 
product A(x)ty[P(x, t)], to be zero in the limits x — > ±oo. 

An important result associated with nonlinear FPE's 
is the H-theorem and its generalizations [H, IqI— Tl3[| . In 
the case of a system subjected to an external potential, 
the H-theorem leads to a well-defined sign for the time 
derivative of the free-energy functional, 



F = U-~fS 



U = 



dx <j>{x)P{x,t) , (2) 



where 7 represents a positive Lagrange multiplier and the 
entropy may be considered in a very general form as, 

f°° d 2 a 
S[P}=J g[P(x,t)]dx; <?(()) = <?(!) = 0; ^|<0 , (3) 

with the condition that g[P(x, t)] should be at least twice 
differentiable. Considering the FPE flTJ, for dF/dt < 



EH ! we obtain, 



d 2 g[P] _ n[P] 
MP] 



dP 2 



(4) 



A relevant outcome of Eq. Q is that the ratio J2[P]/*[P] 
determines an entire class of FPE's associated with a 
single entropic form [t| [l3j . Here we are interested on 
the following type of nonlinear FPE [l4, 15|: 



dP _ d[A{x)P] Dv d_ 
dt dx dx 



P 



_ X 3P 
dx 



(5) 



where D is a constant, v is a real number and the func- 
tionals in Eq. ([1]) correspond to ^[P(x,t)] = P(x,t) and 
Q[P(x,t)] = Dv[P(x,t)] v ~ l . By substituting these quan- 
tities in Eq. Q, integrating and using the conditions (J3|), 
one obtains Tsallis entropy 0, Qjl 17 1, for which, 



g[P] = k 



[P(x,t)] v - P(x,t) 



(6) 



where k = D/~f. Considering an external force A(x) = 
—ax (a > 0), the solution for Eq. ([5]) with initial condi- 
tion P(x,0) = 8(x) is given by the distribution, 



P(x, t) = B(t)[l + /8(t)(l - v)x 2 ] 1 l {v - 1) 



(7) 



where [y\+ =y, for y > and zero otherwise. The time- 
dependent parameters B(t) and f3(t) are defined in such 
a way as to preserve the norm and the form of the distri- 
bution for all times 14. EH. Equation © corresponds 
exactly to the distribution obtained through extremiza- 
tion of Tsallis entropy Eq. (j6|), for the energy defined as 
in Eq. ((2|) , and under standard constraints of probability 
normalization. It is important to notice that, by substi- 
tuting v = 2 — q in Eq. I|T]). one recovers precisely the 
usual distribution of nonextensive statistical mechanics, 
known as g-Gaussian, as obtained through a more general 
definition for the internal energy (l7j |. 
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Additionally, the way in which the system evolves dy- 
namically, as well as its stationary state, provide dis- 
tinctive signatures of anomalous diffusion behavior. For 
small times t <C 1, it is possible to show that the diffu- 
sion propagation front x(t) dictated by Eq. ([5]) advances 
as (x(t)) oc t 2 /(^+i) [:|. In the limit t -> oo, the system 
approaches the stationary state given by, 



P st (x)=B*[l-p* 



„2il/(^-l) 



(8) 



depending on the initial time t , f3* = f3(t Q ){B* / B(t )} 2 , 
and B* = [aB(t ) 2 /2Div/3(t )] 1/(1+,/) . One should stress 
that this form of stationary distribution holds for any 
confining potential <t>(x\ by simply replacing in Eq. ([8]) 
the term x 2 with <f>(x) |14| . 

Next we show that the microscopic behavior of a sys- 
tem of interacting particles undergoing overdamped mo- 
tion is fully compatible with a continuum nonlinear diffu- 
sion equation. Moreover, we find that such a continuum 
formulation for a highly dissipative system corresponds 
to the FPE ([5]) with v = 2. We start by considering 
the equation of motion for a particle i in a system of N 
overdamped particles, 



(9) 



where u, is the velocity of the ith particle, p is the effec- 
tive viscosity of the medium, the first term on the right 
accounts for the interactions among particles, F e (fi) rep- 
resents an external force, and r\ corresponds to an un- 
corrected thermal noise with zero mean and variance 
{i] 2 ) = ksT/fj,. Here we consider a short-range repul- 
sive particle-particle interaction in the form, J(r) = 
G(\f\/X)f, where f is the unit vector along the axis con- 
necting each pair of particles, and A is a characteristic 
length of the short-range pairwise interaction. 

To obtain a continuum description of this system Q , 
we perform a coarse graining of Eq. ([9]) , starting from the 
Fokkcr-Planck equation for the probability distribution 
of the particle coordinates V(fi, fjsr, t), 

^ = E V*(-/> + h B TV{P) , (10) 

i 

where ft is the force on the particle i given by Eq. (|9|). 
By introducing the single particle density p[f, t) = 
(Y]a S 2 (r — f*j)), where the average is made over the dis- 
tribution V(fi, r/v, t), one obtains, 



dp 



dVj(r-r ,/ )p (2) (^' : ".*) 



+F e (r)p(r, t)] + k B TV 2 p{f, t) , (11) 



where p^ 2 \f, r ', t) is the two-point density. If we 
now assume that the approximation p^> (r,r' ,t) ~ 



p(f,t)p(f',t) is valid, we can then coarse grain the 
particle-particle interaction force to obtain Q, 



d 2 r'J{r-r')p(r' ,t)~-aVp{r,t) ; a = d 2 rr-J(f)/2 



(12) 

where only length scales larger than A were considered. 

Here we investigate the motion of particles in a two- 
dimensional narrow channel of size L x x L y under an 
external force in the x-direction, F e = — A(x)x. We then 
assume that the concentration is only weakly dependent 
on the transverse y-coordinate, p{r,t) ~ p{x 7 t). After 
collecting all force terms, we obtain, 



/' 



dp 
di 



0_ 

dx 



dp 
dx 



A(x) 



d 2 P 
dx 2 



(13) 



Interestingly, by introducing the probability P(x, t) — 
(L y /N)p(x,t), defining D = (aN)/(i/L y ), and setting 
/i = 1 and T = 0, we recover precisely the FPE ([5]) for 
v = 2. Herein we restrict our study to a restoring har- 
monic force, A(x) — —ax (a > 0). In the limit T = 0, 
the steady-state solution of Eq. (fT3")) is given by 



psi 



2a 



(x 2 e - x 2 ); \x\ < x e , 



(14) 



with x e = (3Na/2aL y ) 1 / 3 . This result is identical to 
Eq. ([8]) if there one adopts v = 2 and (i* —x~ 2 . 

We now show through direct molecular dynamics (MD) 
simulations that the behavior of a typical overdamped 
system at T = is fully compatible with the solution 
(|14p , therefore representing a microdynamical realization 
of Tsallis thermostatistics. As an example, we consider 
a system of vortex lines moving on a type II supercon- 
ductor substrate 0, , but the results presented here 
could nevertheless be extended to other forms of short- 
range repulsive interactions In this case, the effective 
viscosity is given by p — ^qH^/luc 2 , where $o is the 
magnetic quantum flux, c is the speed of light, uj is the 
resistivity of the normal phase, and H C 2 is the upper criti- 
cal field. Moreover, the commonly adopted vortex- vortex 
interaction is J(r) = /o-Ki(|r|/A)r , where K\ is a mod- 
ified Bessel function decaying exponentially for \r\ > A, 
the pre-factor is given by f = & 2 /(8ir\ 3 ) and A corre- 
sponds to the London penetration length [21| . 

We perform MD simulations with N = 800 flux lines 
placed at random within a two-dimensional channel of 
sizes L x = 100A, L y = 20A, where periodic boundary 
conditions are imposed in the y-direction. In these sim- 
ulations, we first set T = and a confining external 
force is applied in the x-direction with restoring constant 
a = 10 _3 /o/A. The equations of motion (0) are nu- 
merically integrated and the system evolves in time until 
a stationary state of mechanical equilibrium is reached, 
which is identified here in terms of an invariant density 
profile. As shown in Fig. [lj the obtained stationary pro- 
file of particle density is clearly parabolic in shape, in per- 
fect agreement with the theoretical prediction of Eq. (|14l) . 
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FIG. 1: Profile of the density of particles at stationary state 
and T = 0, obtained from molecular dynamics by integrating 
Eq. © (empty circles), as compared to the theoretical esti- 
mate Eq. (| 14[) . The best fit to the data of a quadratic function 
gives the parameter a = 2.41/oA 3 (full line). The position x 
is measured in units of A, whereas the steady-state density 
p(x) is expressed in units A~ 2 . 



The least-squares fit to the simulation data of a quadratic 
function results in the parameter a « 2.41/oA 3 . Strictly 
speaking, the predicted value a = ir fg\ 3 , calculated from 
Eq. (|12[) . is valid only for unconfincd systems, if we as- 
sume that the density p varies slowly within the inter- 
action range of a particle. The discrete character of the 
interacting particles therefore leads to a correction on this 
theoretical prediction. In the absence of external forces, 
i.e., for a = 0, our results show that the front propa- 
gation at early times evolves as (x(t)) oc i 1/3 @. Once 
more, the compatibility between this dynamical scaling 
and the anomalous diffusion behavior intrinsically asso- 
ciated with Eq. ([5]) for v — 2 confirms the validity of our 
approach. 

Having shown that overdamped particles in the limit 
of T — 0, typified as interacting flux lines on a type II su- 
perconductor substrate, obeys Tsallis statistics with an 
cntropic index v = 2, we now analyze the effect of fi- 
nite T on this system. From Eq. (fT3|). one can envisage 
the competition between two types of diffusion, which 
are associated, respectively, with the strength of interac- 
tions between vortexes, controlled by the parameter, a, 
or equivalcntly, D, and the temperature of the thermal 
bath, T. In this way, the ratio k B T/a plays a crucial 
role in the time evolution of the system. As for the case 
T = 0, in the absence of external forces, i.e., for a = 0, 
the diffusion behavior of the system for fcgTCa should 
be governed by the anomalous features associated with 
the index v = 2, namely, (x 2 ) oc t 2 / 3 , whereas normal 
diffusion prevails for fc^T ;g> a, i.e., (x 2 ) oc t. In the 
presence of a restoring external force and for T > 0, a 
stationary-state analytical solution for Eq. (jT3")) can still 



0.5.- 
0.4 '. 
^ 0.3 - 
a. 0.2 . 
0.1 '. 



0.0 
0.5 




T=0.1 



—i — i — i — i — i — i — i — 

O MD 

— Fokker-Planck 

T=2.0 




-200 -100 
X 




100 200-200 -100 
X 



100 200 



FIG. 2: Comparison between the density of particles pro- 
files in the stationary state obtained from molecular dynam- 
ics (empty circles) and the theoretical prediction Eq. (|15|) (full 
lines), for different values of the temperature of the thermal 
bath, T. The position x is measured in units of A, the density 
is in units A -2 , and T is given in units (/oA/fcs). The param- 
eter a = 2.41/oA 3 is the same adopted for T = 0, whereas p(0) 
is determined by conservation of the total number of particles 
in the system, p(x)dx — (N/L y ). For low values of 

T, e.g., T = 0.1/oA/fcB, the density profile is approximately 
parabolic. At intermediate values of T, the system can be de- 
scribed by a W-Lambert function Eq. (|15l) . At high values of 
T, the profile becomes typically Gaussian, as illustrated from 
results at T = 4/ A/fcs. 



be obtained, 
p(x) = 



ap(0) 



ax 



k B T 2k B T 



(15) 



where the M^-Lambert function is defined implicitly 
through the equation W(z)e w ^ = z (see [13] and refer- 
ences therein). In order to test this prediction, extensive 
MD simulations have been performed for different values 
of T. As depicted in Fig. [21 we find very good agree- 
ment between the continuum model solution Eq. (|15|) . 
and MD results of density profiles at steady-state. For 
consistency, we consider the same numerical value used to 
adjust the profile at T — 0, namely, a = 2.41/ A 3 , while 
the parameter p(0) in Eq. (|15p has been determined by a 
conservation constraint in the total number of particles, 
J^ 2 /2 p(x)dx=(N/L y ). 

The effect of increasing T is to gradually change the 
density profile from a parabolic to a Gaussian shape. 
Indeed these two limits correspond, respectively, to the 
particular cases of Tsallis distribution with index v = 2 
and the standard Boltzmann-Gibbs statistics. Next we 
show that the profiles obtained from Eq. (|15p . for finite 
T, are in fact associated with a novel type of entropy 
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form. By comparing Eqs. (fT3|) and (QJ, wc obtain that 
#[P(»,i)] = P(a;,<) and ft[P{x,t)} = 2DP(x,t) + k B T. 
Substituting these quantities into Eq. Q, integrating 
twice and using the conditions ([3]), one gets, 



D 

1 



[P(x,t)-P 2 (x,t)] 



k B T 

7 



[P(x,i)\nP(x,t)] 



(16) 



where 7 is a positive Lagrange parameter defined through 
a free-energy functional like the one given by Eq. (|2|), 
and a; is a conveniently rescaled variable. This functional 
leads to the following entropic form: 



S[P] = -[1- / dxP 2 {x,t)] 
J-00 

/oo 
dx P(x,t)]nP(x,t) . (17) 
-OO 



k B T 



7 

Equation (|T7|) is precisely the sum of Tsallis entropy with 
v = 2, which appears as a consequence of many-body 
interactions, and Boltzmann-Gibbs entropy, which comes 
from the thermal noise of the bath. The entropy is zero 
only for D = T = 0. One may now extremize this mixed 
entropy under the constraints, /_ P(x,t)dx = 1 and 
U = J™ [4>{x) — (f>o]P(x, t)dx, by defining the functional, 



X[P(x,t)) 



7 
k B T 



-a 2 



S[P] 



"1 



dx P(x, t) 



dx [(j){x) - <j) ]P(x,t) . (18) 

The extreme condition ST,/SP = finally leads to 

P(x)exp[^ P(x) = 

exp{-(l + ai)-a 2 [#c)-0o]} , (19) 

which may be written in terms of IF-Lambert's function, 
i.e., in the same form as Eq. (fT5|) . by considering a har- 



monic potential 



; /2, multiplying both sides by 



2D/(k B T), identifying 02 = l/(k B T), defining z as the 
right term in Eq. (JTSJ), andff(z) = 2DP(x)/(k B T). This 
confirms that the entropy Eq. (fTT)) is directly related to 
the stationary solution Eq. (JTSJ of the FPE (|13|). 

In summary, we have shown through a nonlinear 
Fokker- Planck formalism and confirmed through MD 
simulations that a system of interacting particles under- 
going overdamped motion at T = 0, where T is the tem- 
perature of a thermal bath connected to the system, can 
be considered as a physical realization of Tsallis ther- 
mostatistics with an entropic index v = 2 in Eq. ([6]). At 



high values of T, the classical Boltzmann-Gibbs behavior 
is recovered. At intermediate values of T, our approach, 
also confirmed by MD simulations, leads to a stationary 
solution for the corresponding nonlinear FPE that can be 
expressed in terms of the VF-Lambert function. As a con- 
sequence, we disclose a novel mixed entropic form based 
on a linear combination between Tsallis and Boltzmann- 
Gibbs entropies. 
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